//      speed.cpp
//
//      Copyright 2009 Stefan Fruhner <stefan@kandinsky>

//#define DEBUG

#include <iostream>
#include <cstdlib>
//#include <MarchingCubes.h>
#include <vector>

// includes from mylibs
#include <mylibs/mymesh.hpp>
#include <mylibs/mymatrix.hpp>
#include <mylibs/mystring.hpp>
#include <mylibs/mymesh.hpp>
#include <mylibs/lookuptable.hpp>
#include <mylibs/myinifiles.hpp>
#include <mylibs/cmdline.hpp>
#include <mylibs/gipl.h>

#include <myImg/IGBheader.hpp>

using namespace std;
using namespace mylibs;
using namespace mylibs::cmdline;

#ifndef print
#define print(x) (cout << __FILE__ << ":" << __FUNCTION__ <<"() (" <<__LINE__ << ") : " << x << endl)
#endif

#ifndef printe
#define printe(x) (cerr << __FILE__ << ":" << __FUNCTION__ <<"() (" <<__LINE__ << ") "<< "Error : " << x << endl)
#endif

bool verbose	 = false;
mystring igbfile = mystring();
mystring mesh 	 = mystring();
mystring output	 = mystring();
double time_val 	 = 0.;
double resolution	 = 100.;
//double tol		 = 1.e-99;

// Points that define a line
Point p1;
Point p2;

template <typename T>
void process_data(const T *data, ostream &o);

template <typename T>
int read_IGBfile(IGBheader& h, ostream &o){
	int     items = h.x()*h.y()*h.z();
	int slicesize = items * h.data_size();


	int slices = (time_val - h.org_t()) / h.inc_t();		// compute the number of slices to seek
	gzseek((gzFile) h.fileptr(), 1024 + slices*slicesize, SEEK_SET); // seek to the desired time slice, (the size of the header is 1024 bytes)

	double t = h.org_t() + slices * h.inc_t();
	cmdline::msg(" - time value : " + toString(t) + " " + h.unites_t());

	T *data = new T[items]; 								// allocate memory
	int numread = gzread((gzFile) h.fileptr(), data, slicesize); 	// read one slice of data
	if (numread != slicesize) return -1; 					// if no data was read stop here

	// do whatever you like with the data

//	int idx = ..; 											// get the 1D index of the voxel in the mesh
//	T &value = data[idx];							 		// get the value
//	o << i_t <<"\t"<< t << "\t" << value << endl;
	process_data(data,o);
	// ... done
	delete[] data;

	return numread;
}

int main(int argc, char** argv){
//!####################### Load settings from file #############################"

	myIniFiles ini(argc, argv);

	ini.set_info("space_series ... Funktionswerte entlang einer Geraden [p1;p2] interpolieren.");

	ini.register_param("mesh"   , "m", "mesh file in CARP or Tetgen style");
	ini.register_param("output" , "o", "filename of the output (default: <input basename>.dat)");
	ini.register_param("igbfile", "d", "filename of the igb file to load");
	ini.register_param("time"   , "t", "Time value to extract data for. Default: 0.");
	ini.register_param("resolution"   , "r", "Resolution of the data to be written.");

//	ini.register_param("tol",NULL, "tolerance value , usally a small number. Default: 1.E-5");

	ini.register_param("p1", NULL, "Point p1: starting point of the line");
	ini.register_param("p2", NULL, "Point p3: end      point of the line");

	ini.register_flag("verbose", "v", "Additional information on stdout.");
	ini.check();


//	tol 		= ini.read("tol", tol);
	mesh 		= ini.read("mesh"	,	mesh	);
	output 		= ini.read("output", output		);
	igbfile	 	= ini.read("igbfile",	igbfile	);
	time_val	= ini.read("time", 		time_val);
	resolution 	= ini.read("resolution", resolution);

	verbose = ini.exists("verbose");
	vector<double> pt1;	ini.readVector("p1", pt1);
	vector<double> pt2;	ini.readVector("p2", pt2);
	p1.set(pt1);
	p2.set(pt2);

//! ########################## Error handling ##################################

	if (mesh.empty()	) cmdline::exit("Need a mesh to run ... " );
	if ( igbfile.empty()) cmdline::exit("Need an igb file to run ... " );

//! ############################### main program ###############################

	filebuf fb;
	string fn;
	{
		if (verbose) cmdline::section("Reading mesh and igb file");

		gzFile file = gzopen(igbfile.c_str(), "r");
		if (not file) cmdline::exit("Could not read igb file.");

		IGBheader h(file);	// read in igb header

		if (verbose){
			TELL(h.x());
			TELL(h.y());
			TELL(h.z());
			TELL(h.t());

			print(" - x = "<< h.org_x() <<" + i_x * " << h.inc_x());
			print(" - y = "<< h.org_y() <<" + i_y * " << h.inc_y());
			print(" - z = "<< h.org_z() <<" + i_z * " << h.inc_z());
			print(" - t = "<< h.org_t() <<" + i_t * " << h.inc_t());

			print(" - pt1 = " << pt1);
			print(" - pt2 = " << pt2);
		}

		gzseek(file,1024, 0); // seek to the end of the header

		mystring fn(igbfile);
		if (output.empty()){
			if (fn.lower().endswith("igb.gz")) fn = fn(0, -7).file_base(true);
			else								 fn = fn.file_base(true);
			fn.file_new_ext(".dat");
		} else {
			fn = output;
		}

		fb.open(fn.c_str(),ios::out); // opening file to write
		ostream os(&fb);			  // create an ostream for that file

		switch (h.type()){
			case IGB_BYTE     : read_IGBfile<unsigned char >(h, os); break;/* -- byte ------------*/
			case IGB_CHAR     : read_IGBfile<char           >(h, os); break;/* -- Char ------------*/
			case IGB_SHORT    : read_IGBfile<short          >(h, os); break;/* -- short -----------*/
			case IGB_LONG     : read_IGBfile<long           >(h, os); break;/* -- long ------------*/
			case IGB_FLOAT    : read_IGBfile<float          >(h, os); break;/* -- float -----------*/
			case IGB_DOUBLE   : read_IGBfile<double         >(h, os); break;/* -- Double ----------*/
			case IGB_INT      :	read_IGBfile<int             >(h, os); break;/* -- integer ---------*/
			case IGB_UINT     :	read_IGBfile<unsigned int   >(h, os); break;/* -- unsigned integer-*/
			case IGB_COMPLEX  :	/* -- 2 x float (real part, imaginary part) - */
			case IGB_D_COMPLEX:	/* -- 2 x Double (real part, imaginary part)  */
			case IGB_RGBA     :	/* -- 4 x byte (red, green, blue, alpha) ---- */
			case IGB_STRUCTURE:	/* -- Structure ----------------------------- */
			case IGB_POINTER  :	/* -- void * -------------------------------- */
			case IGB_LIST     :	/* -- List   -------------------------------- */
				printe("Datatype not implemented : " << h.type()<< " " << h.TypeID2Name(h.type()));
				break;
			default:
				printe("Datatype unknown : " << h.type() << " " << h.TypeID2Name(h.type()));
				exit(-2);
				break;
		}

		fb.close();

//		mystring xl(h.unites_t());
//		if (xl == "mv") xl = "mV";
//		xl = "time / " + xl;
//		if (xlabel.empty())	xlabel = "time / " + mystring(h.unites_t());
//		if (ylabel.empty())	ylabel = "data / " + mystring(h.unites() );
//		mystring yl(h.unites()); yl = "{/Symbol F}_e / " + yl;
//		write_gnuplot_file(fn);


		gzclose(file);
	}

//! ############################### finalization ###############################
	ini.save();

	return 0;
}

template <typename T>
void process_data(const T *data, ostream &o){

	// load the mesh
	SurfaceMesh m(mesh);

	o << "# " << time_val << " , # elements " << m.elements() << endl;
	o << "# Data values along the line between p1 and p2, with "<< endl;
	o << "#  p1 = " << p1 << endl;
	o << "#  p2 = " << p2 << endl;
	o << "#  resolution  = " << resolution << endl;
	o << "#  original file : " << igbfile << endl;

	Point dir = p2 - p1;	// length of the line

	size_t steps = dir.abs() / resolution;

	dir.normalize();
	Point pt(p1);
	// go along the line
	for (size_t i = 0; i < steps; ++i){
		// interpolate value for the actual point
		try{
			T value = m.interpolate(pt, data);
			// save value to the file
			o << resolution * i << "\t" << value << endl;

		} catch (SurfaceMesh::Exception_InterpolationFailure &e) {
			cerr << e.what() << endl;
		}

		pt += resolution * dir; // select next point
	}
}
